function [tsOut,IT] = detrendZeroMeanMinimal(ts,NC)
%detrendZeroMeanMinimal - Computation of detrended, zero mean time series, 
%minimal approach
%
%       [tsOut,IT] = detrendZeroMeanMinimal(ts,NC)
% 
% This function computes the detrended, zero mean version of the time
% series related to the tsData object ts for the number of components
% NC (allowed values: 2 for East-North, 3 for East-North-Vertical)
% considering a very simple straight line model. 
% If the detrending is successful, the output variable tsOut is upgraded
% with the obtained detrendedZeroMeanE,..., properties and IT is true.
% If the detrending cannot be carried out, tsOut = ts and IT is false.
% If NC is undefined, empty or is neither 2 nor 3, NC = 3 is used.
%
% See also CMEstackFiltering, detrendZeroMean.

% G. Teza, 2022.

if (nargin < 2) || isempty(NC) || (numel(NC) > 1) || ~ismember(NC,[2,3])
    NC = 3;
end

if NC == 2
    VC = [0 1];
else
    VC = [0 1 2];
end

tsOut = ts;
for h = 1:NC
    if VC(h) == 0
        if ~isempty(ts.cleanedE) 
            t = ts.cleanedE.t;
            p = ts.cleanedE.p;   % data in m
            offs = ts.offsetsE;
        else
            IT = false;
            return
        end
    end
    if VC(h) == 1
        t = ts.cleanedN.t;
        p = ts.cleanedN.p;   % data in m
        offs = ts.offsetsN;
    end
    if VC(h) == 2
        if ~isempty(ts.cleanedV)
            t = ts.cleanedV.t;
            p = ts.cleanedV.p;   % data in m
            offs = ts.offsetsV;
        else
            fprintf('WARNING: station %s, vertical data unavailable \n',...
                ts.statName);
            continue
        end
    end
    lts = numel(t);
    tf = serial2frac(t);    % to have t in fractional year
    if isempty(offs)
        [~,~,~,v,~,~,~] = wlinfit(tf,p);     % velocity in m/y
        et = v*(tf-tf(1));  % raw trend
        pd = p-et;          % raw residual
        mpd = mean(pd);
        pdm = pd-mpd;       % zero mean residual time series
    else
        nos = numel(offs);
        pdm = zeros(lts,1);
        for k = 1:(nos+1)
            if k == 1
                t1k = t(1);
            else
                [~,I1k] = min(abs(t-offs(k-1)));
                if t(I1k) <= offs(k-1)  % the offset could be not in t
                    t1k = t(I1k+1);
                else
                    t1k = t(I1k);
                end
            end
            if k == nos+1
                t2k = t(end);
            else
                [~,I2k] = min(abs(t-offs(k)));
                if t(I2k) < offs(k)  % the offset could be not in t
                    t2k = t(I2k);
                else
                    t2k = t(I2k-1);
                end
            end
            dtk = (t1k:t2k)';
            Itk = ismember(t,dtk);  % true for the t in interval dtk
            tfk = tf(Itk);
            pk  = p(Itk);
            if numel(tfk) > 3
                [~,~,~,vk,~,~,~] = wlinfit(tfk,pk);
                etk = vk*(tfk-tfk(1));  % raw trend
                pdk = pk-etk;           % raw residual
                mpdk = mean(pdk);
                pdmk = pdk-mpdk;        % detrended, zero-mean in [t1k t2k]
                pdm(Itk) = pdmk;
            end
        end
    end
    if VC(h) == 0
        tsOut.detrendedZeroMeanE.t = t;
        tsOut.detrendedZeroMeanE.d = pdm;
        tsOut.detrendedZeroMeanE.method = 'Minimal';
    end
    if VC(h) == 1
        tsOut.detrendedZeroMeanN.t = t;
        tsOut.detrendedZeroMeanN.d = pdm;
        tsOut.detrendedZeroMeanN.method = 'Minimal';
    end
    if VC(h) == 2
        tsOut.detrendedZeroMeanV.t = t;
        tsOut.detrendedZeroMeanV.d = pdm;
        tsOut.detrendedZeroMeanV.method = 'Minimal';
    end
end
IT = true;